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ABSTRACT 

We present a generalisation of surface photometry to the higher-order moments of the line-of- 
sight velocity distribution of galaxies observed with integral-field spectrographs. The gener¬ 
alisation follows the approach of surface photometry by determining the best fitting ellipses 
along which the profiles of the moments can be extracted and analysed by means of harmonic 
expansion. The assumption for the odd moments (e.g. mean velocity) is that the profile along 
an ellipse satisfies a simple cosine law. The assumption for the even moments (e.g velocity dis¬ 
persion) is that the profile is constant, as it is used in surface photometry. We test the method 
on a number of model maps and discuss the meaning of the resulting harmonic terms. We 
apply the method to the kinematic moments of an axisymmetric model elliptical galaxy and 
probe the influence of noise on the harmonic terms. We also apply the method to SAURON 
observations of NGC 2549, NGC 2974, NGC 4459 and NGC 4473 where we detect multiple 
CO- and counter-rotating (NGC 2549 and NGC 4473 respectively) components. We And that 
velocity profiles extracted along ellipses of early-type galaxies are well represented by the sim¬ 
ple cosine law (with 2% accuracy), while possible deviations are carried in the fifth harmonic 
term which is sensitive to the existence of multiple kinematic components, and has some anal¬ 
ogy to the shape parameter of photometry. We compare the properties of the kinematic and 
photometric ellipses and And that they are often very similar, but a study on a larger sample 
is necessary. Finally, we offer a characterisation of the main velocity structures based only on 
the kinemetric parameters which can be used to quantify the features in velocity maps. 

Key words: methods; data analysis - techniques: photometric -techniques: spectroscopic - 
galaxies: kinematic and dynamics -galaxies: photometry 


1 INTRODUCTION 

Over the last three decades broad band observations of early- 
type galaxies were successfully analysed by a method commonly 
called surface photometry or, simply, photometry. This method is 
based on the analysis of isophotal shapes of the projected sur¬ 
face brightness. The development of the method was stimulated 
by the empirical discovery that the isophotes of early -type galax¬ 
ies are reproduced by ellipses to better than 1 per ce nt jKenl|l984|: 
[LaMj fTgSSbt IPavis et^IlQS.'j : I.Tedrzeiewskill98^ IPeletieret^ 
1199(1) . Although the isophotes are elliptical in shape to high ac¬ 
curacy, under careful examination, many isophotes of early-type 
galaxies d o show differences from pure ellipses at a level of « 
0.5% (e.g. iBender et aljri988l : IPeletier et alJll990l) . The true suc¬ 
cess of photometry was the ability to measure these deviations 
and clas sify early-type galaxies accordingly into disky and boxy 
objects iUaue 3 |198.5 j : I.Tedrzeiewsldlll98^ (Bender & Moellenhofill 
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Il987h . When combined with information on total luminosity and 
spatially resolved spectroscopy, it followed that the duality of 
photometric properties of early-type galaxies is reflected in a 
duality of kinematic properties, where faint disky objects were 
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). High resolution imag- 


lUauer et ^Il99.5t iRest et all 1200 ll) . again based on photometric 
analysis, revealed new properties of early-type galaxies that deep¬ 
ened the division of the galaxies in two groups. These remark¬ 
able set of discoveries resulted in a re vised classification scheme 
of galaxies iKormendv & BendeJl99(il) . 


It is fair to say that our increased knowledge of early-type 
galaxies (partially) comes from photometry and the ability of the 
method to harvest and describe in a compact way the informa¬ 
tion from two-dimensional images. However, the integrated light 
remains a limited source of information about the internal struc- 
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ture and thus the true nature of galaxies. Kinematic information 
is also essential to fathom the complexity of galaxies and, es- 
pecially, two-dimensional kinematic information is required (e.g 


Franx et alJl^lStatleJl99lHStatler & Fi^l994irtatleill994al 

Two-dimensional velocity maps were, until recently, only 
possible for objects with clear emission lines. Such maps were 
the r esult of stud i es of . e.g., FlI with radio interferometers, 
(e.g. JSwthgrs_etalJ_ | 2002 h . 

(e.g. 


CO with millimetre interferometers 
Ha with Fabry-Perot spectrographs 
te.g. iHemandeTTt al.ll200.5h . The advent of integral-field spectro¬ 


graphs (e.g. TIGER, OASIS, SAURON, PMAS, CMOS, SINFONI, 
OSIRIS), however, has brought two-dimensional kinematic mea¬ 
surements to classical optical and near-infrared wavelengths. Here 
it is possible to probe both the stellar absorption- and gas emission¬ 
lines, which may co-exist in the same potential with very dif¬ 
ferent spatial distributions and dynamical structures. The wealth 
of features seen in st ellar kinematic maps of early-type galaxies 
iEmsellem et alJ2004 confirms the usefulness of two-dimensional 
data, but also poses a problem to efficiently harvest and interpret 
the important features from the maps. 

An approach using harmonic expansion was developed for 
the analysis of two-dimensional velocity maps of disk galaxies. 
This method divides a vel ocity map into in dividual rings (the so- 
called tilted-ring method, [Segem^ ^ 8 ^_ar^_perfonns a_Jia^ 
moni c expansion along t hese rings te.g. [Sinne \1 119781: iTeuberj 
19^ ; iFranx et'^ 1 199^ ISchoenmakers et aiT^gTOwon^^l 
2004h . This analysis, however, is based on the assumption that the 
emission-line emitting material is confined to a thin disk struc¬ 
ture. Clearly, a spheroidal distribution of stars typical of early-type 
galaxies does not have the same dynamical properties as a gas disk. 
To explore those intrinsic properties, one needs a more general 
method which will not be based on assumptions about the nature 
of the observed system (e.g. a disk), but will rely solely on the 
properties of the investigated observables (e.g. surface brightness 
or velocity). 

Here we explore such a general method. Photometry earned 
its spurs in intensive applications over the last three decades and 
offers a natural starting point for the analysis of two-dimensional 
kinematic maps, although one can, in principle, invent a number 
of basis functions which can describe a two-dimensional distribu¬ 
tion. In this paper, we present a generalisation of photometry to 
the higher-order moments of the line-of-sight velocity distribution 
(LOSVD). This generalisation is based on the theoretical fact that 
surface brightness is itself a moment of the LOSVD, and on our em¬ 
pirical discovery that the velocity maps of many early-type galax¬ 
ies can be well reproduced by a simple cosine law along sampling 
ellipses. We call our method kinemetry^, which reduces to photom¬ 
etry for the investigation of surface brightness distributions. 

In Sectionl^we present the theoretical background and moti¬ 
vation for the new method. Section^presents the technical aspects 
of the method. The meaning of the kinemetric coefficients and their 
diagnostic merits for different model maps are presented in Sec- 
tion|3 In Section0we apply the method on kinematic maps of an 
axisymmetric model galaxy, present the application of kinemetry to 
actual observations and characterise typical structures on velocity 
maps. We summarise our conclusions in Section|^ 


2 THEORETICAL BACKGROUND AND MOTIVATION 


The dynamics of a collisionless stellar system is fully specified by 
its phase-space density or distribution function f = fix, v, t) (e.g. 
iBinnev & Tremainall987^ . but this quantity is not measurable di¬ 
rectly. When observing external galaxies, we measure properties 
that are integrated along the line-of-sight (LOS). An additional 
complication is that the galaxies are viewed from a certain direc¬ 
tion, and we actually observe only those projected properties of the 
integrated distribution function. 

The observables, which reveal only averages over a large 
number of unresolved stars, are the surface brightness and the 
LO SVD, sometime s also called the velocity profile^ (for a review 
see Ide Zeeuwll994i . The surface brightness is given by 

ki-{x,y)=[ dz [ [ [dv fif,v). (1) 

7 los J J J 

The rest of the observable information is carried by the LOSVD, 
which relates to the distribution function via 


£.iv;x,y)= [ dz [ [dv,^dvy f {f, v), (2) 

7 los j j 

where (x,y,z) are the three spatial coordinates, oriented such that the 
LOS is along the z-axis. From these equations it is trivial to show 
that surface brightness is the zero-th order moment of the LOSVD. 
Higher-order moments, such as (U), or (U^), can easily be related 
to the observables of the LOSVD: mean velocity V, velocity dis¬ 
persion a and higher-order moments, commonly parametrised by 
Gaus s-Hermite coefficients ivan der Marel & Fran^ 19^ iGerhaid 
^9^, hi and (14 being the most used. 

The kinematic moments of stationary triaxial systems show 
a high degree of symmetry which can be expressed through their 
parity. The mean velocity is an odd moment, while the velocity 
dispersion is an even moment. In practice, this means that a two- 
dimensional map of a given moment shows corresponding symme¬ 
try. Maps of even moments are point-symmetric, while maps of odd 
moments are point-anti-symmetric. In polar coordinates on the sky 
plane this gives: 


pL^ir,e + 'K) = Peir,e), 

p,oir,9 -\-It) = -fioir,9), (3) 


where pe and po are arbitrary even and odd moments of the 
LOSVD, respectively. Furthermore, if the observed system is ax¬ 
isymmetric, the even moment of the LOSVD will also be mirror- 
symmetric or, correspondingly, an odd moment will be mirror-anti- 
symmetric: 


Peir,Tt-9) = Peir,9), 

Poir,n-9) = - Pair, 9). (4) 

A kinematic moment that satisfy both symmetries (eqs.|3and|3 is 
said to be bi-(anti)-symmetric. The existence of these symmetries 
can be used to simplify the harmonic analysis. Point-symmetric 
moments give rise only to even harmonic terms, while point-anti- 
symmetric moments will require only odd harmonic terms. Mirror- 
(anti)-symmetry additionally requires that there is no change in the 
moment’s position angle. For a more detailed discussion on the in¬ 
fluence of these equations on the terms of the harmonic expansion 
see Appendix I aI 


^ This name was introduced by|Co 2 in^ '^l200i who presented a prelim¬ 
inary discussion on this topic. 


^ In this work, however, the term ‘velocity profile’ (or generally, kinematic 
profile) will be used to denote a set of velocity (kinematic) measurements 
extracted from a map along a certain cuiwe. 
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Following these conventions, surface brightness, as the zeroth 
order moment, is an even moment. As the surface brightness and the 
kinematic moments of the LOSVD are moments of the same distri¬ 
bution function, it is natural to analyse them similarly, and we can 
ask in which way one can generalise the method of photometry to 
work on the higher moments of the LOSVD. Considering symme¬ 
try, the even moments do not require much change in the method, 
but the odd moments need a new working assumption. Also, while 
the isophotes of integrated light are elliptical to a high accuracy, 
the choice of sampling curves is not as obvious for odd kinematic 
moments. 

Since the light distribution of early-type galaxies is elliptical 
it is reasonable to assume that an expansion along ellipses would 
also be suitable for the extraction of profiles from kinematic mo¬ 
ment maps giving some insight in their structure and properties. 
The choice of the actual ellipticity, however, is somewhat arbitrary. 
For example, it is possible to expand along the galaxy isophotes 
or along ellipses that correspond to deprojected circles. In the first 
case, one follows the photometry and probes regions of the same 
projected surface brightness. In the second case, one samples lo¬ 
cations with equal intrinsic radii, but it is necessary to assume an 
inclination for the galaxy, which is usually difficult to estimate. It 
may also not be physically justified to use a constant ellipticity for 
the whole velocity map, as the ellipticity of the light distribution 
may change with radius. 

As a special case, it is possible to expand a kinematic 
map a l ong circ l es. Th is ap proach was i nvesti gated bv ICopin et al.l 
i200ll) . ICopinI i2002h and lKrainovi3 i2004l . and although it is 
straight-forward and requires no a priori assumption, a large num¬ 
ber of harmonic terms are often necessary, preventing a simple 
description of the maps and complicating the analysis (see Ap¬ 
pendix]^ for a discussion). Additionally, an expansion along el¬ 
lipses can r educe to expansion along circles in limiting cases (see 
section l4~3l . lFranx et in <1989l) also note the advantages of extrac¬ 
tion along circles, even for photometry, but conclude that an ap¬ 
proach based on ellipses is superior. 

Consider the first moment of the LOSVD, the mean velocity. 
In specific cases, the geometry of the system offers a choice for 
curves along which velocity profiles can be extracted. If a velocity 
map traces material in a thin disk, such as is often the case in HI, 
CO and Ha gas emission maps, velocity can be well represented by 
circular motion. If the inclination of the disk is known, a circular 
orbit in the plane of the galaxy can be projected as an ellipse on the 
sky. A velocity profile extracted along this ellipse is then of simple 
cosine form, which can be expressed by 

y(i?, ■!/)) = Vb + Vc(f?) sinicost/), (5) 

where R is the radius of a circular ring in the plane of the galaxy 
(or the semi-major axis length of the ellipse on the sky), Vb is the 
systemic velocity, Vc is the ring circular velocity, i is the ring in¬ 
clination (i = 0 for a ring seen face-on) and tp is the azimuthal 
angle measured from the projected major axis in the plane of the 
galaxy. The inclination i is related to the axial ratio or flattening, 
q, of the ring’s major and minor axes. For a given flattening, the 
ring’s velocity profile is most similar to the velocity of a circular 
orbit with radius R and inclination cos i = q. In other words, the 
flattening defines an ellipse on the sky, with ellipticity e = 1 — q, 
which corresponds to a circle describing the orbits in the plane of 
the galaxy. 

If one assumes that the velocity profiles extracted from veloc¬ 
ity maps (either stellar or gaseous) of early-type galaxies can be 
described by the simple cosine law of eq. 01°^ high accuracy, 


then the first step in the analysis of velocity maps is to determine 
those ellipses along which the velocity profile is best described by 
eq. In analogy with photometry, the next step would be to de¬ 
termine the deviations of the velocity profile from the cosine law, 
which can be achieved by the higher-order Fourier analysis. 

We discovered that the velocity profiles extracted from the 
maps of the SAURON sample can indeed be well described by the 
cosine law of eq. Note that this observation suggests that the 
velocity maps of early type galaxies closely resemble the observed 
kinematics of inclined circular disks. One has, however, to take a 
step back and recognise that spheroids are not disks, and the two 
cases should be clearly differentiated during the interpretation of 
the results. The fact that the velocity maps of galaxies are well ap¬ 
proximated by the thin disk model allows us to develop a simple 
generalisation of photometry based on eq. 0, but also reflects an 
interesting property of the internal structure of galaxies, which is 
worth exploring further (see section 15^ for a few examples while 
the results for the rest of the SAURON sample will be reported in a 
future paper). 

The method presented in this paper is applicable to all mo¬ 
ments of the LOSVD. The quality of kinematic data, however, usu¬ 
ally decreases rapidly with increasing order of the kinematic mo¬ 
ment: the signal is generally simply too weak for a good enough 
extraction of the higher-order moments from raw data. Perhaps it is 
even fair to say that only the mean velocity maps reach the quality 
of photometric observations from thirty or so years ago, when sur¬ 
face photometry was introduced. For these reasons, as well as for 
presentation purposes, we restrict our analysis to the mean velocity 
maps. The method, however, can be straight-forwardly applied to 
the higher order moments (see Section lsTTI for an example). 

3 THE METHOD 

In this section we present the details of the new algorithm extending 
surface photometry to the maps of odd moments of the LOSVD. 

3.1 Harmonic expansion 

Fourier analysis is the most straight-forward approach to charac¬ 
terise any periodic phenomenon. A velocity map K (a, ip) can be 
divided into a number of elliptical rings yielding velocity profiles 
which can be described by a finite (and small) number (N + 1) of 
harmonic terms (frequencies): 

N 

K{a,ip) = Ao(a) -F ^ A„{a) 3 m{nip) + Bn{a) cos{nip), (6) 

n=l 

where ip is the eccentric anomaly, that in case of disks corresponds 
to the azimuthal angle of eq. a is the length of the semi¬ 

major axis of the elli ptical ring. 

As presented bv ICarteJ il978h . lK^ il983h and ijedrzeiewskil 
^9^), to find the best sampling ellipse along which a profile pe 
should be extracted, one minimises the harmonic terms sensitive to 
the parameters of the ellipse (centre, position angle and ellipticity 
for a given semi-major axis length). In the case of an even moment 
(e.g. surface brightness) the best fitting ellipse parameters are given 
by minimising Ai, Bi, A2 and B2 coefficients of the harmonic ex¬ 
pansion. This means that a photometric ellipse is determined when 
a sufficiently good fit to pe(a, ip) = Ao{a) is obtained. In the case 
of an odd moment of the LOSVD (e.g. mean velocity), given our 
empirical finding, the sampling ellipse parameters are determined 
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by requiring that the profile along the ellipse is well described by 
eq. 0, or, generally, /ro(a, t/;) = Bi{a) cos ip. 

Ellipse parameters (position angle, centre and flattening) can 
be determined by a minimising a small number of harmonic terms. 
In Appendix I bI we investigate which terms are sensitive to changes 
of the ellipse parameters. We find that: 

- an incorrect kinematic centre produces non-zero Aq, A 2 , B 2 , 
and to a lower level Ai, A 3 and B 3 coefficients, 

- an incorrect flattening (inclination) produces a non-zero B 3 
coefficient, 

- an incorrect position angle produces non-zero Ai, A 3 and B 3 
coefficients. 

In the study of disk velocity maps. l^hoenmakers et 
in their Appendix A2, discuss the same sensitivity, but for small 
deviations. In this l imit, our results reduce to the analytically de¬ 
rived conclusions of ISchoenmakers et alJ^1997^ . It is, however, im¬ 
portant to note, that whatever the combination of incorrect centre, 
flattening or position angle, only the n < 3 harmonic terms are 
affected. 


3.2 The algorithm 

Given our discovery that the cosine law along sampling ellipses de¬ 
scribes well velocity maps of galaxies, a small number of harmonic 
terms is needed to determine the best fitting ellipse. In this case 
eq. reduces to 

K{ip) = Aq + Al sm{ip) + Bi cos('i/)) 

-I-A 2 sin(2')/)) -I- B 2 cos{2ip) 

+A3 sin(3t/)) -I- B3 cos(3t/)) (7) 

where ip is the eccentric anomaly. In practice, as it is done in pho¬ 
tometry, we sample the points uniformly in ip and the ellipse coor¬ 
dinates are then given by Xe = o cos ip and ye = qa sin ip. This 
implies that, if the ellipse is projected onto a circle, the sampling 
points are equidistant. The selection of K(ip) values is obtained by 
a bilinear interpolation of the o bserved map, also in the same way 
it is usually done in photometry iledrzeie^ 

The parameters of the best sampling ellipse for an odd kine¬ 
matic moment are obtained by minimising: 

= Al + Al + Bl+Al + Bl (8) 

from eq. 0. This non-linear fit is performed in two stages. The 
first stage consists of minimising on an input grid of flattenings 
q and position angles F. The second stage repeats the fit using a 
non-linear least-square curve fitting algorithm, but now also fitting 
for the coordinates of the centre, taking as initial values q = qmin 
and r = Fruin from the first stage and an estimate of the centre 
(xo,yo). In this way, the true x^ minimum can be robustly deter¬ 
mined and, hence, the parameters of the best sampling ellipse: cen¬ 
tre (Xc,yc), flattening (ellipt icity) and position angle. We use the 
MINPACK implementation iMore et alJll98lll of the Levenberg- 
Marquardt method®. This two stage approach was implemented to 
improve the robustness of the method, preventing the minimisation 
program from choosing possible secondary minima. 

The choice for the length of semi-major axis of consecutive 
ellipses is set by the decrease of surface brightness of the observed 

® We used an IDL version of the code written by Craig B. Markwardt and 
available from: http://astrog.physics.wisc.edu/~craigm/idl 


object. The kinematic maps are necessarily spatially binne d (e.g. 
using adaptive Voronoi binning iCaPDellari & CoDiril2003h . as in 
the examples of Section |3, and as the bins typically increase in 
size with radius, the distance between the rings has to be increased 
to avoid correlation between the ring parameters. The correlation 
should also be avoided in the centre where the bins are usually iden¬ 
tical to the original pixels. For these reasons, we adopted a combina¬ 
tion of linear and logarithmic increase of the semi-major axis length 
given by this expression: a = n -\- l.l'*, where n = 1, 2, 3.... 

The final ellipse parameters obtained by minimisation are then 
used to describe an elliptical ring from which a kinematic profile is 
extracted and expanded onto the harmonic series of eq. where 
the coefficients {An, Bn) are determined by a least-squares fit with 
a basis {1, cos ip, sin ip, ..., cos Nip, sin Nip}. 

The harmonic series can be presented in a more compact way, 

N 

K(a,ip) = Ao{a) + ^ kn{a) cos[n{ip - (pn{a))]. (9) 

n=l 

where the amplitude and phase coefficients (fc„, cpn) are easily cal¬ 
culated from the An, Bn coefficients: 

kn = A‘i + Bl and cpn = arctan ■ (10) 

The sin and cos terms have different properties and, in prin¬ 
ciple, describe different characteristics of the maps. The prop¬ 
erties of the corresponding coefficients have been studied and 
described for the case of thin gaseous disk s ipranx et*^ 1 19941 : 
ISchoenmakers et ^Il997t I Wong et"^ l2004|) . In a general case 
however, namely for the analysis of kinematic maps of triaxial (stel¬ 
lar) systems, those properties do not apply (e.g. the bi term rep¬ 
resents circular rotation only in the thin disk approximation). For 
these reasons, we find it more instructive in this paper to combine 
the coefficients of the same order and characterise the global prop¬ 
erties of kinematic maps. In what follows, although the implemen¬ 
tation of the method determines individual coefficients (and can be 
straight-forwardly used for the special cases of disks), they will be 
presented in the spirit of eq. 

3.3 Internal error estimates 

The uncertainties of the resulting parameters can be estimated from 
the measurement errors of the input kinematic data. The ellipse 
parameters determined by the Levenberg-Marquardt least-squares 
minimisation code have formal one-sigma uncertainties computed 
from the covariance matrix. If an input parameter is held fixed, or 
if it touches a boundary, then the error is reported as zero. Simi¬ 
larly, the harmonic terms obtained in the linear least-squares fit in 
the second part of the procedure have their formal one-sigma errors 
estimated from the diagonal elements of th e corresponding co vari- 
ance matrix, as described in Section 15.6 of IPress et al.Ul992h . 

We performed Monte Carlo simulations to estimate the ro- 
bustness of t he ob tained internal errors. We used the SAURON 
jBacot^^^lJj20^^ velocity map of NGC 2974 presented in 
lEmseiie^^ta^^2004^ and its measurement errors to construct one 
hundred realisations of the velocity map. The velocity in each bin 
was taken from a Gaussian distribution of mean equal to the ob¬ 
served velocity and standard deviation given by the corresponding 
measurement error. Typical errors on the velocities of the SAURON 
observations are about 5 km s“ ®. All realisations provide a distribu¬ 
tion of values from which we estimate one-sigma errors. The Monte 
Carlo uncertainties are in a good agreement with the formal errors. 
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Figure 1. Four model maps representing different but typical velocity maps of triaxial eai'ly-type galaxies. Velocity maps were constructed of two components 
(see text for details). The position angle of the first (nuclear) component changes from left to right: 0, 45, 90 and 180°. The size of bins increases towards the 
edge of the maps, as expected form observations, to keep constant signal-to-noise ratio. Note also the different opening angles of the two components on the 
maps. 


Note, however, that these errors are statistical errors only and do not 
include systematic effects. We investigate the resulting accuracy of 
the kinemetric coefficients in Section|3 


4 KINEMATIC PARAMETERS AND THEIR MEANING 

The crucial aspect of any analysis method is the interpretation of 
its results. In this section, we describe in more detail the way the 
method works and demonstrate the meaning of the resulting kine¬ 
metric parameters. 

4.1 Model velocity maps 

We constructed a set of model velocity map s combining two c om- 
ponents using the circular velocity from the lHemauislI i 199(1) po¬ 
tential, weighted with two different surface brightness distributions. 
The Hemquist potential was chosen because it approximates rea¬ 
sonably well the density of real early-type galaxies. Its circular ve¬ 
locity is given by 14 = VGMrf{r + ro), where G is the grav¬ 
itational constant, M is total mass of the system, and ro is scale 
length of the potential. The VGM factors for the two components 
were 850 and 1500 (in units of kms“^ arcsec^/^), respectively, 
while ro scale lengths were set to 5”and 15”, defining a central and 
large scale component, respectively. 

The light distribution of the first component was described by 
an exponential law typical of disks, and that of the second com¬ 
ponent by an r^^^ law typical of spheroids. They were both nor¬ 
malised to a maximum value of unity. The light distribution scale 
lengths, re, of the two components were set to 7”and 15”, re¬ 
spectively. The components were inclined (60° and 40°, respec¬ 
tively) and combined with 4:1 relative weights, respectively, ro¬ 
tating the central component for different position angles T = 
(0°, 45°, 90°, 180°) in order to construct different features on ve¬ 
locity maps possible to occur in triaxial early-type galaxies. 

The final result are four model velocity maps, which satisfy 
the symmetry relations of eq. presented in Fig.Q Each model 
has a distinct nuclear component with a somewhat higher maxi¬ 
mum velocity and a smaller opening angle than the large scale, 
spheroidal component. We define the opening angle loosely as a 
quantity that describes the tightness of the iso-velocity contours. It 
is discussed in more detail below. The main difference between the 
models is in the orientation of the nuclear component. All model 
maps were constructed to simulate realistic observations with the 
SAURON integral-field spectrograph, where the size of the spatial 


bins increases from the centre outwards to maintain constant signal- 
to-noise ratio. 

4.2 Kinemetric analysis of model maps 

The parity of an odd kinematic moment guarantees that the even 
terms in the expansion vanish in the absence of measurement errors 
and assuming that the centre is given correctly. These conditions are 
satisfied for the model velocity maps and, while applying the algo¬ 
rithm of section lT^ we simplified eqs. Q and by not minimis¬ 
ing the contribution of A 2 and B 2 . However, during the harmonic 
analysis we expanded the extracted velocity profiles including two 
additional odd harmonic terms, A 5 and B 5 , which describe the de¬ 
viations from circular rotation (see eq.|5}- 

A special case is that of the zeroth-order term given by the 
Aq coefficient. It is not shown in the subsequent analysis because 
the absolute level of the maps was set to zero by construction. For 
velocity maps, Ao gives the systemic velocity of each ring. In the 
case of even kinematic moments, the meaning of Ao is different: it 
is the dominant term describing the amplitude of the moment (e.g. 
intensity for surface brightness, amplitude for velocity dispersion 
of velocity dispersion maps, see section IsTl . 

Figure]^ shows the kinemetric parameters which describe the 
model velocity maps. All parameters are plotted versus the semi¬ 
major axis length of the ellipses. The first and second row show the 
orientation (F) and flattening (q) of the ellipse along which the ve¬ 
locity profile was extracted. These parameters uniquely specify the 
best fitting ellipse, along which the contribution of Ai, A 3 and B 3 
terms is minimal. The following rows show plots of the harmonic 
coefficients given in the compact form of eq. < 10 > : ki and fcs/fci. 

The meaning of these parameters can be associated with the 
visible structures on the maps as follows: 

(i) r is the kinematic position angle. Like the photometric po¬ 
sition angle, it describes the orientation of the velocity map, and 
it is related to the projection of the angular momentum vector (see 
Appendixlcl. F traces the position of the maximum velocity on the 
map. 

(ii) q is the axial ratio or flattening of the ellipse and it describes 
the flattening of the map. It can be interpreted as a measure of the 
opening angle of the iso-velocity contours. Maps that have large 
opening angles appear round (large values of g), while maps with 
small opening angle appear (small values of q). In the special 
case of a disk, where the motion is confined to a plane (spiral galax¬ 
ies, gaseous disks, etc), the flattening is directly related to the incli- 
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model 0 model 45 model 90 model 180 



Figure 2. Kinemetric coefficients of the model maps. From left to right the plots show the results for model 0, model 45, model 90 and model 180. From top 
to bottom: kinematic position angle F, flattening q, coefficients of harmonic expansion fci and k^jki. Note that the sudden rise of fcs/fei in model 180 at 9”is 
due to the nearly zero value of the rotation velocity fci at that radius. 


nation of the disk, q = cos(i), where the projected ellipses corre¬ 
spond to intrinsically circular orhits. 

(iii) fci is the dominant coefficient on maps of odd kinematic 
moments. For velocity maps, it describes the amplitude of bulk mo¬ 
tions (rotation curve). It contains contributions from both cos tp and 
sin Ip terms, and in the special case of a thin disk geometry, the Ai 
and Bi terms should be considered separately. In the case of perfect 
circular motions (Ai and other higher terms are zero), fci describes 
the full motion on the map. In other words, if higher order coeffi¬ 
cients are zero (or negligible), then the velocity map is identical to 
the velocity map of a disk viewed at inclination i. 

(iv) fcs is the combination of harmonic coefficients that are min¬ 
imised in the fitting process. This coefficient is sensitive to the, cen¬ 
tre, flattening and position angle of the ellipse and in special cases it 
will be larger then zero indicating a velocity profile with a complex 
structure (multiple extremes) originating from the overlap of mul¬ 
tiple kinematic components with different orientations (see Fig.|3- 
This coefficient, however, does not bring any additional informa¬ 
tion to q and F and we do not show it. On the other hand, if q and 
F are fixed and not fitted, as is often the case in studies of disk ve¬ 


locity maps, then this term brings the first correction to the simple 
rotational motion. 

(v) fcs is the first higher-order term which does not define the 
parameters of the best fitting ellipse. It represents the higher-order 
deviations from simple rotation and points to complex structures on 
the maps. This term is thus sensitive to the existence of separate 
kinematic components. It is a kinemetric analogous of the photo¬ 
metric term that describes the deviation of the isophote shape from 
an ellipse (diskiness and boxiness) 

Comparison of the structures on the model maps (Fig.0 and 
the radial dependence of the kinemetric coefficients (Fig.|3 illus¬ 
trates the descriptive power of the method. Our model maps were 
built of two components which had different flattening and kine¬ 
matic position angle. It is the role of the kinemetric coefficients to 
recognise these characteristics and describe them (Fig.|3- 

The components of model 0 were aligned, but the nuclear com¬ 
ponent iso-velocities had a small opening angle, while those of the 
large-scale component had a large opening angle. This is reflected 
in the different flattenings q. In the inner 10”, the disky component 
dominates and the map is rather while beyond this region the 
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Figure 3. Velocity profile taken from the model 90 map along a circle 
(q = 1) at t=9" , in the transition region between the two components. Over¬ 
plotted with diamonds is the circular velocity profile at the same position. 
The observed profile exhibits a complex structure which can not be repre¬ 
sented by requiring that the third-order terms in the harmonic expansion be 
zero. The local extremes at t/i Ri 100° and 280° originate in the nuclear 
component, while the global maximum and minimum are given by the large 
scale component. 


spheroidal component dominates and the map looks more round. 
In the regions where only one component dominates, the recovered 
axial ratio corresponds correctly to the inclination of that compo¬ 
nent. While the fci terms show the amplitude of the rotation, the 
kz and fcs coefficients are close to zero. This is expected because 
the models were built of two maps of essentially circular velocity. 
A small kz terms exists in the transition region between the com¬ 
ponents. This suggests that the combination of aligned but separate 
velocity components will give rise to non-zero kz. 

A confirmation of this effect is given by model 180, where the 
two components are aligned but counter-rotating. In this configu¬ 
ration the components are clearly visible on the q and ki panels 
and kz rises strongly in the region coinciding with the change of 
rotational sense (change in the kinematic position angle by 180°). 

Differences in kinematic position angle between the compo¬ 
nents, as in the maps of model 45 and model 90, generate more 
power in higher order coefficients. The position angle is recovered 
correctly in the regions dominated by one component, while in the 
region between the components it traces the position of (combined) 
maximum velocity. It is the same with the flattening, which follows 
the change of the map topology, reaching extreme round values 
in the transition region where both components equally dominate. 
The reason for this is the complex shape of the velocity profile, for 
which it becomes impossible to make the Az and Bz terms along 
any ellipse equal to zero (Fig.|3- In such a case, the fitting routine 
chooses the equal ratio of axes as this is the most general case, con¬ 
firming the robustness of the approach. The two stage fit, direct grid 
fitting followed by non-linear minimisation, is necessary to provide 
the required robustness in such regions. 

The extent of the transition region depends on the relative posi¬ 


tion angle, larger angles resulting in a larger transition region. This 
is reflected in the behaviour of the higher-order terms, which be¬ 
come more important and influence a larger area with increasing 
misalignment of the components. The increasing misalignment is 
effectively separating the components, creating a complex, multi- 
peaked velocity profile such as in Fig.|3 Note that strong contribu¬ 
tion of higher-order terms in this examples is not necessarily repre¬ 
sentative of real galaxies, as the features on the velocity maps were 
not modelled to represent specific observed cases. 

Figurel^shows the two-dimensional representation of the anal¬ 
ysis. For each model, we reconstructed a map using only the bi term 
which is not minimised, and we named those new maps ’circular 
velocity’. We also reconstructed a map using all 6 terms of the har¬ 
monic expansion, which we named ’kinemetric velocity’. The dif¬ 
ference between the two velocity maps shows the contribution of 
the higher order terms {kz and kz). These residual maps of model 0 
and model 180 show distinct five-fold symmetries coming from the 
high kz term. Residuals of other models show a combination of the 
kz and kz terms in the transition region where velocity profiles are 
similar to the one in Fig.|3 

The two examples with five-fold symmetries are interesting 
when compared with the corresponding model maps. Both resid¬ 
ual maps have positive amplitudes on the negative side of the ma¬ 
jor axis. This is, however, not the case for the model maps. The 
model 0 map has two components which both rotate in the same 
direction, while for model 180 the inner component counter-rotates 
with respect to the outer component. Both models have a strong kz 
term, and hence five-fold residual structure in the transition region. 
Interestingly, the amplitude of the residuals of model 180 at e.g. 
r = 7" along the major axis is negative while the mean rotation at 
the same radius is positive, belonging to the region dominated by 
the counter-rotating component. It is thus possible to infer that the 
residuals actually trace the contribution of the main body rotation in 
this region, although this is not visible on the velocity map. A simi¬ 
lar situation is visible in model 0, but the two components co-rotate 
and the residuals show the same sense of rotation. 

We conclude this exercise stating that kz is sensitive to the 
existence of different components, where the increase of kz in the 
presented models is correlated to larger misalignments between the 
components. 

4.3 Special cases 

There are two possible cases of velocity maps for which the de¬ 
termination of kinemetric coefficients will be degenerate. The first 
case arises from the effec t of solid body rotation. As shown in 
ISchoenmakers et aOjim, while it is possible to measure the po¬ 
sition angle r(a), it is impossible to uniquely determine the flatten¬ 
ing q{a), the centre and the systemic velocity. This degeneracy is 
present because all iso-velocity contours are parallel, and the open¬ 
ing angle is 180°. 

Another extreme case occurs when rotation is very slow. Such 
velocity maps are often found in giant ellipticals which are thought 
to be triaxial bodies. Here the velocity profiles are roughly constant 
and do not show the parity expected from an odd moment of the 
LOSVD. There is no gradient in velocity and as a consequence the 
flattening and position angle cannot be determined. 

In these two cases the new method cannot be applied blindly 
because it is not possible to interpret harmonic terms in the way 
we showed up to now. There are two possible ways to proceed. One 
way is to make an assumption about the flattening and the kinematic 
position angle (e.g. use photometric flattening), and continue with 


© 2005 RAS, MNRAS OOO. mrTTl 





8 


Davor Krajnovic et al. 


V model 



-ZO -1[} 0 10 20 


— 150 km/s 150 


V circular 




-20 -10 0 10 20 


-150 km/s 150 


V klnemetry 





-20 -10 0 10 20 


-150 km/s 150 


V residual 



-20 -10 0 10 20 
anew; 


— 15 km/5 15 


Figure 4. Model velocity maps and the results of the kinemetric analysis. From top to bottom: maps comesponding to model 0, model 45, model 90 and model 
180 (as noted at the end of each row). From left to right: model maps, maps of circular velocity, reconstructed model maps using 6 harmonic terms, residual 
between circular and reconstructed velocity maps. Overplotted on the model maps are the best fitting ellipses along which the velocity profiles were extracted 
and analysed. Note the change in flattening and position angle of the ellipses. 


the harmonic analysis of the velocity profiles extracted along those 
ellipses. 


A more general way, for both cases, is to extract the velocity 
profiles along concentric circles and perform the harmonic analy¬ 
sis on those. In this way one does not bring any assumption in the 
analysis, but it is still possible to determine the rotation velocity and 
kinematic position angle (if present) and detect trends in higher- 
order harmonic terms. A detailed discussion of extraction along cir¬ 
cles and the meaning of the coefficients is given in AnnendixlXI An 
automatic transition towards circles can be implemented by adding 
a penalty term to the to bias the solution towards round shapes, 
when the gradient is small (or iso-velocities are parallel). The good¬ 
ness of fit is then determined by Xp = + OiP{q). We also 

note that an expansion on circles might be good for the analysis of 
higher-order moments of the LOSVD which normally have rather 
poor signal-to-noise ratio. 


5 APPLICATION 

In this section we apply our kinemetry formalism to three different 
examples. Our wish is to justify the choice of expansion on ellipses 
along which the profiles satisfy the simple cosine law. We first apply 
the method to the odd and even kinematic moments of an axisym- 
metric model of a galaxy. Then we apply the method to a few ob¬ 
served galaxies. Finally, we offer a characterisation of structures on 
the velocity maps which can be recognised through the behaviour 
of the kinemetric parameters. 


5.1 Axisymmetric two-integral model 

The working assumption of the method introduced in Section 
is that there exists an ellipse along which the profile of a kine¬ 
matic moment can be described by a simple cosine law (eq. 1^. 
This pro file is easily associated with the circular motions in a thin 
disk (e.g. ISchoenmakers et which show regular velocity 
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Figure 5. Maps of mean velocity (right), velocity dispersion (middle) cr, and Gauss-Hermite moment /13 (left) for the model galaxy with added intrinsic scatter 
corresponding to a realistic measurement with SAURON. The overplotted ellipses on the velocity map are the surface brightness contours. The plots in each 
right and left column represent the kinemetric parameters (F, q, fci, fcs/fci) describing the maps. The plots in the middle column are kinemetric analogs of 
the standard photometric parameters: position angle, flattening, intensity (Aq term) and the shape parameter (S 4 ). Horizontal line in the last plot of the middle 
column marks the division between boxiness and diskiness of iso-cr contours. Blue lines in the middle and right column plots correspond to the noiseless maps 
of (T and hi, whose contours are overplotted on the cr and hi map. 


maps. More general stellar systems with axial symmetry, however, 
also show very regular velocity maps. Here we demonstrate that 
our working assumption of eq. is satisfied in realistic models for 
axisymmetric galaxies. 

We constructed an axisymmetric model galaxy with a distri¬ 
bution function which depends only on the energy and one com¬ 
ponent of the angular moment um. The model was constructed us- 
ing the [ Hunter & (^)ianl il99.3h contour integration as in method 
lEmsellen ^^iiy^^^) to resemble the SAURON obser vations of 
NGC 2974 iEmsellem et al 12004 l&ai no vie et all200.5h . The dis¬ 
tribution function of the model yields the full LOSVD from which 
the observable kinematics can be extracted. This was done on a 
large 100" x 100" field by fitting a Gauss-Hermite series (V, a, 
hi and hi) to spatially binned data resembling typical SAURON ob¬ 
servations. In order to mimic real observations, we also assigned 


typical SAURON observational errors to each model observable (the 
relevant errors here are on velocity, velocity dispersion and hi, re¬ 
spectively 5 kms“^, 7 kms“^ and 0.03). These errors were used 
to add intrinsic scatter to the noiseless model data by means of a 
Monte Carlo realisation. F or further details of the model, we refer 
the reader to Section 5.1 of lKrainovic et in] 11^. Model maps of 
velocity, velocity dispersion and hi are presented in Fig.0 

We performed a kinemetric analysis of two odd and one even 
moment of the model galaxy. The left column of Fig.|5|presents the 
results for the velocity map. The flattening is constant over most 
of the map with a small but continuous change in the centre. The 
kinematic position angle is very well determined and remains con¬ 
stant. The first harmonic term, fci, describes the rotational motion 
in the model galaxy and peaks around 15". The higher order term 
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ks, mostly consistent with zero, shows some deviations on a 0.5-1 % 
level in the centre and slightly more at larger radii. 

The middle column of Fig. [^presents the analysis of the first 
even moment of the LOSVD: the velocity dispersion map. As men¬ 
tioned earlier, because even moments have the same parity as the 
surface brightness (zeroth moment), the kinemetric analysis of an 
even moment, such as velocity dispersion, is analog to the photo¬ 
metric analysis of the surface brightness. We demonstrate this by 
fitting the velocity dispersion map using our code, but with the re¬ 
quirement that = Ao{a) along a trial ellipse. This yields 

the usual photometeric parameters: position angle, ellipticity, inten¬ 
sity and the shape parameters describing the deviations from the el¬ 
lipse. In our case of an even kinematic moment we recognise them 
as the: kinematic position angle T, flattening (q = 1 — e), Aq and 
f?4 harmonic terms, while other harmonic terms in eq. ^ are zero. 
The zeroth harmonic term, Aq, is the dominant term in the expan¬ 
sion and it represents velocity dispersion profile with radius. The 
fourth harmonic term, B4, can be used to describe the deviations of 
the iso-cr contours from an ellipse, where, as in photometry, nega¬ 
tive values are indication of boxiness, while positive values of disk- 
iness. The given velocity dispersion map is in general boxy, but be¬ 
tween 10"and 20”it becomes disky. The same region is marked by 
a change in flattening (elipticity) and orientation (position angle). 

Observed velocity dispersion maps are usually noisier than the 
velocity maps. We analysed a noiseless model velocity dispersion 
map (without intrinsic scatter) to determine with what confidence 
one could analyse an observed velocity dispersion map. The results 
are also plotted in the middle column of Fig.|3 and they show that 
the main features of kinemetric coefficients are well recovered, sug¬ 
gesting that the method could be used to analyse velocity dispersion 
maps obtained from data with high signal-to-noise ratio. 

The right column of Fig.l^presents the results of a similar ex¬ 
ercise, but this time on the observed third moment of the LOSVD 
(given by Gauss-Hermite parameter / 13 ). The noise on this map con¬ 
siderably distorts the structure of the map. The kinematic position 
angle is again very well determined, being, as expected, offset 180° 
from the kinematic angle of the velocity map, and the fci term traces 
well the amplitude of the (13 with radius. The flattening changes 
with a jump at ~ 25”. The higher-order coefficients are however 
substantially influenced by the noise. The fcs term is often more 
than 10 % of fci. We repeated the analysis on the noiseless (no in¬ 
trinsic scatter) map of /13 moment. The new values of q, PA and 
fci, also shown in Fig.|3 are in a good agreement with the values 
extracted from the map with the intrinsic scatter. This is, however, 
not the case for the last term, fcs, where the noiseless model term 
is small (consistent with zero) everywhere except in the region of 
the change of the flattening q. This feature is completely absent 
from the noise dominated fcs term of the intrinsic scatter map. We 
conclude from this exercise that, although the method can be ap¬ 
plied straight-forwardly to the higher moments of the LOSVD, the 
signal-to-noise ratio in current state of the art measurements is gen¬ 
erally not high enough to reliably recover kinemetric parameters of 
the kinematic moments beyond velocity dispersion. 

The negligible higher-order harmonic terms for both odd mo¬ 
ments confirm that the assumption of a simple cosine law for the 
kinematic moments is justified in the case of an axisymmetric 
galaxy. In other words, the method assumes that to zeroth order 
galaxies are axisymmetric rotators. The non-zero higher-order har¬ 
monic terms suggest the presence of multiple components (fcs 7^ 0). 


A change in q and F, however, hints to departures from axisymme- 
try**. 


5.2 SAURON velocity maps 

We now apply t he method to a few galaxies observed in the 
SAURON survey ide Zeeuw et alJl2002h . We selected the velocity 
maps of NGC 2549, NGC 2974, NGC 4459 and NGC 4473 because 
they show diverse velo city structures. The maps were presented in 
lEmsellem et alJ yOO^ and are reproduced in the first row of Fig.|^ 
Inspection of the two-dimensional maps shows that the iso-velocity 
contours of NGC 2549 behave differently in the central and outer 
regions, suggesting the existence of two velocity components. The 
velocity map of NGC 2974 is remarkably regular, with a constant 
opening angle of the iso-velocities. NGC 4459 has a round velocity 
map with a large opening angle, and the velocity map of NGC 4473 
shows an unusual decrease of the velocity along the major axis. 

Our purpose here is not to undertake a detailed study of the 
properties of these galaxies, but to illustrate quantitatively the appli¬ 
cation of the new method. Kinemetric profiles for the velocity maps 
of the four galaxies are shown in Fig. Q We also show the pho¬ 
tometric position angle and flattening (where qphot = 1 — ^phot) 
for a tentative comparison between the surface brightness (as an 
even moment of the LOSVD) and velocity (an odd moment of the 
LOSVD). A detailed study of the photometric properties of the 
SAURON galaxies will be presented in a paper in preparation by 
the SAURON team. We summarise the main results as follows: 
NGC 2974 is largel y co nsistent with axisymm etry (but see 
lEmsellem et alJ ^2003^ and iKrainovic et alJ j2005h for bar signa¬ 
tures), and the simplicity of the velocity map is reflected in the 
kinemetric coefficients: the fcs term is below 1-2%, while the flat¬ 
tening and position angle are both constant over the map except in 
the central few arcsecs where the departures from the axisymmetry 
are the strongest. The photometric and kinematic position angles 
and, especially, the flattenings agree very well. 

NGC 4459 is a kinematically round map with large opening angles 
reflected in the high values of q. The position angle of the veloc¬ 
ity map changes in the central region but remains constant beyond 
7”. In the same region, the velocity reaches a maximum and then 
decreases, with a hint of a rise towards the edge of the map. Sim¬ 
ilarly to NGC 2974, the fcs coefficient is small and below 2% in 
the investigated region. It is possible that there is a separate central 
velocity component, indicated by the change in F and fci. The sim¬ 
ilarity with NGC 2974 is also reflected in the comparison with the 
photometric parameters: the position angles agree very well beyond 
5” while flattenings agree over most of the map. 

NGC 2549 has a multi-component velocity map. Its central region 
(r < 5”) is described by a linear rise in the velocity and a low 
flattening, rising up to 7” — 8” where the velocity reaches its max¬ 
imum. Further away the velocity drops and beyond 10” so does 
the flattening (changing again at r > 17”). These features are also 
recognisable on the velocity map. The behaviour of the fes term is, 
however, not obvious from the map. It is on average larger than 
in the two previous examples over the whole map, and it shows a 
clear rise between 10 and 15”, marking the change between the two 
components. The kinematic and photometric position angles differ 
inside 7” and outside 15”, while flattening of the contours is equal 


^ Strictly speaking, in axisymmetric galaxies F has to be constant and equal 
to the photometric position angle. 
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Figure 6. Velocity maps, reconstructed kinemetric velocity maps and the corresponding residual map (constructed as in Fig.|4} of NGC 2974, NGC 4459, 
NGC 2549 and NGC 4473 (from left to right). The best fitting ellipses are overplotted on the velocity maps. Maps were oriented such that the major axis is 
horizontal, and the north-east orientation of the maps is shown with the arrows. Residuals of NGC 2974 and NGC 4459 are smaller then the five-fold residuals of 
NGC 2549 and NGC 4473. Along the major axis, the residuals of NGC 2549 have the same sense of rotation as the velocity map, while residuals of NGC 4473 
show the opposite sense of rotation. 


only beyond 15”, but a general similarity between the two flatten¬ 
ings is noticeable. 

NGC 4473 is our final example. The decrease in the velocity along 
the major axis coincides with the rise of the flattening and k^. The 
kinematic flattening does not match well the flattening of the photo¬ 
metric ellipses, although the position angles agree very well in the 
region of the velocity decrease. 

Fig-Elshows the reconstructed (kinemetric velocity) and resid¬ 
ual (constructed as in Fig.|3 maps in the second and third row, 
respectively. All velocity maps are well reconstructed with a small 
number of terms. NGC 2974 and NGC 4459 have small residu¬ 
als, and the features on their maps are well described with only ki 
terms. Both NGC 2549 and NGC 4473, however, have a strong fcs 
term. In section l4^ using model maps, we showed that large fcs is 
likely to indicate the existence of two kinematic components. The 
relative orientation of the components can be any, but we draw the 
attention of the reader to the co- and counter-rotating cases pre¬ 
sented by model 0 and model 180. Precisely the same behaviour is 
visible in Figs.|3and0for NGC 2549 and NGC 4473, respectively. 

The two co-rotating velocity components are clearly visible in 
the map of NGC 2549: a central fast rotating component with high 
flattening and an outer slower rotating component with somewhat 
smaller flattening. On Fig. |5|it can be seen that along the major axis 
the contribution of the residuals is in the same direction as the bulk 
rotation, confirming the co-rotation of the components. 

NGC 4473 is not as clear a case. The velocity map (Fig. 
does not show any clear signature of the two kinematic compo¬ 
nents which are inferred from the higher-order kinemetric terms. 
The residual map (Fig.|^, however, shows that the residuals are op¬ 
posite to the bulk rotation, suggesting that the two components are 


actually counter-rotating. This configuration can naturally explain 
the observed decrease of rotation with radius. It is thus likely that 
NGC 4473 is made up of two components: a dominant one respon¬ 
sible for the bulk rotation of the stars, and a lighter counter-rotating 
one (perhaps a disk because of its confined contribution along the 
major axis) which bring s more weight to large rad i i. A co nfirmation 
of this result is given bv ICaPDellari & McDermij i200.5l) . who con¬ 
structed Schwarzschild dynamical models of this galaxy and recov¬ 
ered two counter-rotating components in the distribution function. 

The velocity maps presented in this section are by no means 
meant to be representative of all possible velocity structures. It is, 
however, striking to see how the assumption of the method (eq.|3 
is satisfied to a level of 2 % or better for a number of real galax¬ 
ies. On the other hand, some galaxies do show strong deviations in 
the higher-order terms, corresponding to substructure and multiple 
kinematic components. The correspondence between photometric 
and kinemetric parameters is also interesting. While it is expected 
that the two position angles should agree if a galaxy is axisymmet- 
ric, it is not expected that the flattening of the isophotes and the 
kinematic ellipses will be the same. This is satisfied in the case of 
thin disks, where qphot = Qkin = cos(i), but this relation does not 
necessarily apply to spheroidal systems. The correspondence be¬ 
tween qphot and qkin for spheroidal galaxies is likely related to the 
degree of anisotropy of the galaxies. It may also indicate that these 
galaxies contain stellar disks. 


5.3 Characterisation of features on stellar velocity maps 

It is possible to describe the v ariety of structures pre sent in veloc¬ 
ity maps of galaxies, as seen in lEmsellem et alJi2004h . or observed 
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Figure 7. Application of the method to SAURON observations of NGC 2974, NGC 4459, NGC 2549 and NGC 4473. Each panel shows kinemetric coefficients 
for the velocity maps: the kinematic position angle, axial ratio of the best fitting ellipses, and harmonic terms: fci and k^. Dashed horizontal lines are global 
kinematic position angle determined as the median of F (blue line) and with the method outlined in Appendix^ They are almost always indistinguishable. For 
comparison, the photometric values for the position angle and the flattening of the best fitting photometric ellipses are overplotted (diamonds). 


with other integral-field units, using the results of our analysis. The 
features can be present uniformly over a whole map or can be dif¬ 
ferent between the inner and outer parts of the map. Here we offer a 
characterisation of features on the velocity maps (or the maps them¬ 
selves) based only on the kinemetric parameters. Building on the 
results from the previous sections, velocity features can be sorted in 
at least the following four descriptive groups: 

Disk-like Rotation (DR). This group consists of the simplest and 
featureless velocity maps. A DR map shows clear rotation where 
the amplitude of k\ is substantial, but not necessary constant or ris¬ 
ing. The fcs coefficient is consistent with zero, while the kinematic 
position angle T and the flattening q remain constant. 

Multiple Components (MC). Similar to our model 0, MC maps 
exhibit features which can be related to kinematically different, but 
aligned components. The components are kinematically separate if 
the flattening changes in the transition region, or if it is different for 
each component. Additionally, MC are detected by a rise of the fcs 
coefficients in the transition region and possibly but not necessarily 
by a drop of fci. 

Kinematic Twists (KT). Twists in the iso-velocity contours are 
characterised by a change in the kinematic position angle T. Ad¬ 
ditionally, several characteristics of the MC group can be present 
as well, e.g. fcs > 0 or drop in fci, indicating different kinematic 
components. 


Low-level Rotation (LR). A rather special feature on the veloc¬ 
ity maps occurs when there is no detectable rotation. This effect 
is dependent on the instrumental resolution (rotation is too low to 
be measured), but its origin is in the formation scenario and in¬ 
trinsic distribution of orbits. This group can be characterised by a 
requirement that fci remains below a certain limiting value of fci. 
Maps belonging to this class are characterised by high amplitude of 
fcs, often bigger then fci. Due to the indeterminacy of q and T for 
the LR type, the analysis should be performed along circles, which 
changes the meaning of the harmonic terms. 

It is important to note that the exact details, such as the limiting 
value of fci, or the tolerated change in T, which serve to characterise 
the structures, depend on the instrumental settings and the quality of 
the observed data. However, it is clear that the kinemetric parame¬ 
ters could be used to describe and classify the properties of velocity 
maps. 


6 SUMMARY 

We have presented a generalisation of surface photometry to the 
higher-order moments of the LOSVD, observed with integral-field 
spectrograph. We call our method kinemetry. For even moments of 
the LOSVD, kinemetry reduces to photometry. For odd moments. 
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kinemetry is based on the assumption that it is possible to define an 
ellipse such that the kinematic profile extracted along the ellipse can 
be well described by a simple cosine law. This assumption is satis¬ 
fied in the case of simple axisymmetric rotators. Kinematic profiles 
that deviate from axisymmetry or contain multiple kinematic com¬ 
ponents will be more complex, and the residuals can be measured 
through a harmonic analysis of the profiles. We also find that the 
assumption is well satisfied for a number of observed galaxies. 

Since the profiles are extracted along the best fitting ellipse, the 
number of harmonic parameters necessary to describe the profile is 
usually small. A total of six harmonic terms is generally necessary 
to determine the best fitting ellipse (Aq, Ai, A 2 , B 2 , A 3 and B 3 ). 
The deviations from the simple cosine law are carried in higher har¬ 
monic terms. It is not expected that many more higher-order terms 
will be necessary in the general case, and we stop the analysis at the 
As and B 5 harmonics. 

The most important parameters of the analysis are the kine¬ 
matic position angle T, flattening q of the ellipses, and the harmonic 
terms: fci and fcs, where we combined the terms of the same order. 
We briefly summarise their properties: 

- the kinematic position angle traces the maximum velocity and 
describes the orientation of the map. 

- the flattening is related to the opening angle of the iso¬ 
velocities and describes the different components of the maps, 
which can appear flat or round. In the special case of a thin disk 
it gives the inclination of the galaxy; 

- ki is the dominant harmonic coefficient and describes the am¬ 
plitude of bulk motions. It is related to the circular velocity. Higher 
order coefficients, which we normalise to ki, describe the devia¬ 
tions from simple circular motion; 

- fcs is the first higher-order term which is not fitted and quan¬ 
tifies higher-order deviations from rotational motions. It indicates 
complexity on the maps and is sensitive to the existence of separate 
kinematic components. A photometric analog of this term is the co¬ 
efficient that describes the diskiness/boxiness of the isophotes. 

This method can be straight-forwardly applied to all moments 
of the LOS VD observed from different objects (e.g. early- and late- 
type galaxies). However, the analysis and the intepretation of the 
harmonic terms has to be related to the physical nature of the ob¬ 
served objects. An example are velocity maps of thin disks. In this 
case, the kinemetric flattening and position angle are directly re¬ 
lated to the inclination and orientation of the disk, and the results 
of kinemetry are equal to results of the tilted-ring method. If q 
an d r are kept fixed, kineme try r educes to th e harm onic analysis 
of ISchoenmakers et alJ ^1997^ and IWong et alJ j 20 ^. In this case, 
since the velocity profiles are not extracted along the best fitting el¬ 
lipses, all harmonic terms will carry information about departures 
from the simple circular motion. 

We applied the method on one model and four observed galax¬ 
ies. The model galaxy was a two-integral model of an observed 
axisymetric galaxy. We constructed maps of observable odd kine¬ 
matic moments of the LOSVD: mean velocity, velocity dispersion 
and Gauss-Hermite coefficient /is. We analysed both maps to show 
that the method works well on an axisymmetric galaxy (all higher- 
order terms were small or consistent with zero). We also introduced 
an intrinsic scatter and currently realistic measurement uncertain¬ 
ties to show the difficulties expected analysing the higher-order odd 
moments (e.g. / 13 ) obtained from state-of-the-art integral-field ob¬ 
servations. The position angle, flattening and amplitude of the mo¬ 
ment can be realistically recovered, but the signal-to-noise ratio is 
too low for expansion to higher-order harmonics. 


The four galaxies were taken from the SAURON sample. We 
used them to show the descriptive power of the new method: de¬ 
tection of CO- or counter-rotating subcomponents well hidden in the 
kinematics, which could be detected previously only through de¬ 
tailed modelling. We also showed that there are early-type galaxies 
with rotation velocity described by a simple cosine law. The kine¬ 
metric position angle and flattening of these galaxies are in a good 
agreement with the photometric position angle and flattening. 

Based on the kinemetric parameters we characterise several 
features on the velocity maps: Disk-like Rotation, Multiple Compo¬ 
nents, Kinematic Twists and Low-level Rotation. Other, more spe¬ 
cific and instrument dependant features on the maps can be con¬ 
structed from these basic groups. 

We wish to stress that the method works well on velocity maps 
because many of them satisfy the basic assumption of the method to 
the level of a few per cent. This empirical fact is related to the inter¬ 
nal structure of the galaxies. We developed the method of kinemetry 
to harvest the information from the observed maps of moments of 
the LOSVD probing the nature of galaxies. 

An IDL implementation of the described algorithm is avail¬ 
able on the web address: http://www-astro.physics.ox.ac.uk/~dxk/ 
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APPENDIX A: EXTRACTION ALONG CIRCLES 

The simplest case of kinemetric expansion is by extracting velocity 
profiles along concentric circles. This approach has several advan¬ 
tages over the more general case of ’ellipse fitting’. They are: 

- there is no a priori assumption about the structure of the map 

- extraction of harmonic coefficients is a purely linear problem 
(no degeneracy between parameters) 

- reconstruction of maps is straight-forward even in the most 
complicated cases 

These points, however, are weighted against the fact that the 
number of the extracted harmonic terms which describe a more 
complicated map is large, and often the coefficients cannot be 
clearly associated with the physical properties of the objects. In this 
Appendix, for com plet eness, we briefly summarise the analysis of 
ICnnin et ^ i 200 lli and iKrainovid i2004h . 

Requiring that axial ratio g = 1, kinemetry reduces to Fourier 
expansion along concentric circles via eq. where the coeffi¬ 
cients An, Bn are determined by a least-squares fit with a basis 
{1, cos 9, sin 6 ,cos N9, sin N6}, where 6 is now the polar an¬ 
gle, to which the eccentric anomaly, ip, reduces in this case. De¬ 
pending whether the fit is to an odd or an even kinematic moment, 
the terms in the expansion can be selected such that n is odd or 
even, respectively. In a more compact way, the harmonic series is 
presented by the eq. j^, where the amplitude and phase coefficients 
(kn, (pn) are easily calculated from the An, Bn coefficients follow¬ 
ing eq. fTol . 

The following description of the harmonic coefficients is sim¬ 
ilar to the description of the coefficients obtained by ellipse fitting, 
but differs in some crucial details. 

The zeroth-order term, Ao, measures the mean level of the 
map. For the first kinematic moment, this is equivalent to the sys¬ 
temic velocity of the galaxy. For /13 maps, which have mean level 
equal to zero by construction, Aq will, like other even terms, be 
consistent with zero. 

The coefficient fci gives the general shape and amplitude of 
the odd moment map and is the dominant term. The first-order cor¬ 
rection is given by the next odd term, k^. This term can be named 
the morphology term because it describes most of the additional ge¬ 
ometry of the map. Often, the first two odd terms are enough to 
describe the velocity map, although the kinemetric expansion may 
require higher terms as well, fcs and even k^. 

Connected to the amplitude terms are the corresponding phase 
terms, (pn- They determine the orientation of the map, but their con¬ 
tribution depends on the relative strength of the amplitude terms. 
The first phase coefficient, (pi, gives the mean position angle of the 
velocity map (measured from the horizontal axis, 9 = 0). The an¬ 
gle that (pi measures is the position angle of the ki term. This is, 
in general, slightly dilferent from the positions of the maxima on 
the map, which are also influenced by the contribution from higher- 
order terms. However, it does give the global orientation of the map. 
This phase corresponds to previously defined kinematic position an¬ 
gle F. 

The angle (p 3 is the phase angle of the third harmonic: the next 
significant term. For a small amplitude of fes, its contribution to 
the overall orientation will be small, and the position of the maxi¬ 
mum velocity will be given accurately by 9max = (pi - The phases 
of higher-order terms are interesting because, for an axisymmetric 
galaxy, they satisfy the relation: 

77 TT 

(pi-(pi = — (Al) 
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where n G Z and is the i-th order term in the expansion. 

This condition can be easily derived considering that for an 
axisymmetric velocity map the position of the zero velocity curve, 
the curve along which the velocity is zero on the map, is orthogonal 
to the kinematic angle given by (j)i- This means that K(a, 9) = 0, 
for 6 — (j)i -\- -K12. Neglecting the higher order terms, (ki > fca), 
and substituting 6 = (/)i + 7r/2 into eq. ^ one obtains the result of 
eq. lAH . 

Alternatively, deviations from axisymmetry can be quantified 
following eqs. and <A1> . If K(a, 9 = 4>\ -\- 7 r/ 2 ) = AV 7 ^ 0, 
then one finds: 

- cj)3), (A2) 

Ki Kl 

where we express the relation as a ratio of the dominant term in 
the expansion. This relation quantifies the contribution of the fcs 
term due to departures from axisymmetry and can be generalised 
to other higher terms. In cases of large misalignments between the 
kinematic and photometric position angles, should be replaced 
by the photometric PA in the above equation. 

In many ways, analysis of the even moments is analogous to 
that for the odd moments, taking the proper parity into account. 

The dominant term of the even maps is Aq, and it describes the 
absolute level of the map as a function of radius. The next important 
term in the expansion is fc 2 , and it is the morphology term of the 
even moment maps, which describes features such as elongation 
along the minor axis or “bow-tie” shapes, ofte n seen in observed 
velocity dispersion maps iEmsellem et all2004h . The orientation of 
the morphologically distinct features are determined by the phase 
angles (j >2 of the k 2 term. 

The parity of kinematic moments, expressed by eqs. 0 and 
in section]^ has the following consequences on the harmonic 
expansion along circles: 

- point-symmetry requires only even terms in the harmonic ex¬ 
pansion 

- point-anti-symmetry requires only odd terms in the harmonic 
expansion 

- additional mirror-(anti)-symmetry requires only cosine terms 
in the harmonic expansion 

It follows that kinemetry can easily be used for constructing 
two-dimensional kinematic maps of a given symmetry. This is done 
by fixing certain harmonic terms (odd, even or sine terms) to zero; 
e.g. a bi-anti-symmetric map is created by keeping only odd cosine 
terms in the harmonic expansion. Additionally, if the number of 
harmonic terms in the expansion is small, the reconstructed map 
will be smoothed. In this way, kinemetry can be used for removing 
the higher-order harmonics from the data and effectively filtering 
the map. 


APPENDIX B; INFLUENCE OF INCORRECT CENTRE, 
FLATTENING AND ORIENTATION 

Here we test the sensitivity of harmonic expansion of a velocity 
profile given by eq. Q with respect to a chosen incorrect flattening 
(inclination in a thin disk case), position angle and the centre of the 
ellipse along which the profile was extracted. We constructed a sim¬ 
ple velocity map using the circular velocity form the Hernquist po¬ 
tential and an exponential light profile corresponding to the central 
(disk) component of the model 0. The velocity map had a constant 
flattening of go = 0.5 and kinematic position angle of To = 0°. 



-10 
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Figure Bl. The effect of an incorrect flattening on the coefficients of har¬ 
monic expansion of a velocity profile. Upper panels shows the model ve¬ 
locity map having constant orientation (kinematic position angle) and flat¬ 
tening. The solid ellipse has correct flattening and orientation, while the 
dashed ellipse has incorrect axial ratio and correct orientation. The lower 
panel shows residuals between the extracted velocity profile and the cir¬ 
cular velocity B\ cos('i/;) (solid symbols). Overplotted are contributions of 
A\ sin('0) (red line) and cos(3'0) (blue curve). 
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Figure B2. Same as Fig. IB II where the dashed ellipse has correct flatten¬ 
ing but incorrect orientation. Next to the residuals on the lower panel over¬ 
plotted are curves indicating the contribution of different terms: Ai sin(' 0 ) 
(red curve), 2 I 3 sin(3'0) + S 3 cos(3'0) (blue curve), and a combination 
Ai sin('0) + As sin(3‘i/^) + S 3 cos(3'0) (green curve). 
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Figure B3. Same as Fig.|BjJwhere the dashed ellipse of correct flattening 
and position angle was shifted by 1”horizontally and l”vertically from the 
centre. The lower panel presents the residual and conh'ibutions from differ¬ 
ent terms: Ag (black line), j4i sini/i (red line), A 2 sin{2'i/j) + B 2 cos(2t/!) 
(orange curve) and A 3 sin(3'i/>) -I- B 3 cos(3'0) (blue curve). The combined 
contribution of all these terms is overplotted with the green curve. 


Systemic velocity was set to zero. We first extracted a velocity pro¬ 
file from this map at a = 10” along a pre-designed ellipse. The 
ellipse was first chosen to have an incorrect flattening q' = 0.6, 
but correct kinematic position angle To (Fig. lBU . We repeated the 
exercise selecting the ellipse with incorrect T' = 10° and correct 
q (Fig. lB2> . Each of the profiles were then harmonically expanded 
with eq. o . Finally, we shifted an ellipse, defined by correct values 
of q and F, horizontally and vertically for l''(Fig. lB^ . The corre¬ 
sponding velocity profile was expanded using first N = 6 (even 
and odd) terms in eq. 

Harmonic expansion of a profile extracted along an ellipse de¬ 
fined by the correct centre, flattening and orientation would (by 
construction) yields all harmonic coefficients zero expect the Bi, 
which gives the circular velocity of the map. Hence, on the lower 
panels of Figs. lBlIlBTI and lBSI we plot the difference (residuals) of 
the extracted velocity profile and the Bi cos('i/i) term. 

On Fig. IBII overplotted are the contributions of terms: 
Al sin('i/i) and B 3 cos(3(/i). The A 3 sin(3'i/i) is nearly zero, like the 
ai, and is not shown. The contribution of only the B 3 cos(3t/i) term 
can explain the distribution of residuals, confirming that only B 3 
coefficient is affected when choosing an incorrect flattening. 

On Fig. IB2I overnlotted are the following terms: Ai sin(t/i), 
A 3 sin(3'i/i) -I-H3 cos(3'i/i), B 3 cos(3t/>), as well as the combination 
Al sin(t/i) + A3sin(3'i/i) + B 3 cos(3t/i). All of these terms con¬ 
tribute significantly to describe the distribution of residuals. The 
dominant coefficient is Ai, while A 3 and B 3 contribute in this ex¬ 
ample with about 45% and 15% in a mplitude, respectively. In the 
limit of a small incorrect position angle. l^hoenmakers et 
showed that only Ai and A 3 are affected, but our example shows 
that for larger incorrect position angles the B 3 term has to be con¬ 
sidered as well. 

Fig. IB3I presents contributions for Ag, Ai, A 2 , B 2 , A 3 and 


B3 harmonic terms. To the first order, only Ag, A2 and B2 terms 
are sensitive to i ncorrect centring. The same result was also analyt¬ 
ically derived bv ISchoenmakers et alJ ^9^), in the limit of small 
miscentering. However, for large deviation from the true centre co¬ 
ordinates, contribution from the A3 and B3 terms becomes impor¬ 
tant, while AI remains much smaller. 

We conclude this exercise by stating that the ellipse parameters 
for an odd kinematic moment map can determined by harmonic 
expansion of the Ag, Ai, A2, B2, A3 and B3 terms only. 


APPENDIX C: DETERMINATION OF GLOBAL 
KINEMATIC POSITION ANGLE 

The apparent angular momentum Lp of a galaxy, as a projection of 
the intrinsic angular momentum to the plane of sky, is defined as 
<Franxll988h : 

Lp = J" Tp X Svi-c/fp (Cl) 

where rp is a projection to the plane of sky of the vector r at which 
the mean velocity of stars, v, projects to the observed radial velocity 
vector Vr. E is surface density and the integration is over the whole 
projection plane. As noted bv iFranx et al.l il99lt) . it follows from 
eq. <C1> that the apparent angular momentum is fully specified by 
the observed surface brightness and the observed two dimensional 
velocity map. In addition to that, if figure rotation is absent, the ap¬ 
parent angular momentum is parallel to the projection of the intrin¬ 
sic angular momentum. This means that the projection of the full 
three-dimensional velocity space is reduced to the projection of the 
angular momentum. Reversing the argument, measuring the (lumi¬ 
nosity weighted) kinematic position angle of the two-dimensional 
projection of velocities, it is possible to determine the position of 
the apparent angular momentum. 

The position angle of the apparent angular momentum Tkin 
can be determined from eq. ICH using integral-field observations 
of galaxies. This is the most physical way of determination of 
but it is somewhat dependant on the geometry of the maps®, like the 
relative orientation of the map with respect to the orientation of the 
apparent angular momentum. 

A more robust way of determining F^^i^ is based on a direct de¬ 
termination of the kinematic position angle F. Kinemetry measures 
F as a function of radius, but we are interested in a global (indepen¬ 
dent of radius) value of the kinematic position angle, Fg, because 
of a direct correspondence between Fg and Fkin- Fg can be deter¬ 
mined in a robust way using the information in the whole velocity 
map. For any chosen Fg, we construct a bi-anti-symmetric velocity 
map V'(x,y) with the a:-axis along the given Fg - This is done by 
replacing the measured mean velocity V{x,y) inside each Voronoi 
bin, of centroidal coordinates {x, y), with the weighted average of 
the corresponding velocity in the four quadrants, using linear inter¬ 
polation when needed to estimate the velocity: 

^ y{x,y) + V{x,-y)-Vj-x, y)-V{-x, -y) 

(C2) 

^ We do not discuss here the influence of the size of the maps. Clearly, only 
r kin of the observed part of the galaxy can be determined, because eq. iCll 
is defined over the whole projection plane 
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Figure Cl. Comparison of observed velocity map (left) and bi-anti- 
symmetric map (right) at the best fitting F q using the method outlined in 
the text. We used SAURON observations of NGC 4473. North is up and east 
to the left. 

When any of the four symmetric coordinates was outside the con¬ 
vex hull defined by the bins centroids, we used only the remaining 
values to determined V'. Fig. 1C II shows an example of a bi-anti- 
symmetric map. The true F c was defined as the angle which min¬ 
imises 



(C3) 


In this way it is possible to assign simple error estimates to the 
obtained Fg as the range of angles for which < 9, which 

corresponds to the 3cr confidence level for one parameter. 

We used velocity maps from section 1.5.21 to compute the 
global kinematic position angle. Obtained values for NGC 2974, 
NGC 4459, NGC 2549 and NGC 4473 are 43.2, 280.6, 1.0, 92.5, 
respectively. The results are overplotted on Fig. Q where they can 
be compared with the kinemetric position angle, F, which changes 
with radius describing specific features on the maps. They are in 
a very good agreement as can be seen comparing with the median 
values of the kinemetric position angle. 
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